!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! calculate function coefficients for given nozzle wall function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! wall function is y=aw+bw*x+cw*x^2+dw*x^3
! slope of wall is dy/dx=tan(theta)=bw+2*cw*x+3*dw*x^2
! slope of the slope of the wall is d^2 y/dx^2==2*cw+6*dw*x
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   dda    slope of slope of wall at the attachment point
!   aw,bw,cw,dw    wall equation constant
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   all angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine Wall3c()
   use VariableDef
   implicit none
   real*8::dda
   
!   attachment point
   xa=rtd*sin(thetaa)
   ya=yt+rtd*(1-cos(thetaa))
   dda=xa**2/((rtd**2-xa**2)**1.5)+1/((rtd**2-xa**2)**0.5)
!   constant
   aw=(dda*xa**4-4*da*xa**3-2*de*xa**3-4*dda*xe*xa**3+12*da*xe*xa**2+3*dda*xa**2*xe**2&
      &+6*ya*xa**2-6*da*xa*xe**2-12*xa*xe*ya+6*ya*xe**2)/(6*(xe-xa)**2)
   bw=(de*xa**2+dda*xe*xa**2-2*da*xa*xe+da*xe**2-dda*xa*xe**2)/((xe-xa)**2)
   cw=(2*da*xa-2*de*xa-dda*xa**2+dda*xe**2)/(2*(xe-xa)**2)
   dw=(de-da+dda*xa-dda*xe)/(3*(xe-xa)**2)
   
   ye=aw+bw*xe+cw*xe**2+dw*xe**3   ! exit lip point radical location

end subroutine Wall3c
